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ABSTRACT 

The known extrasolar multiple- planet systems share a surprising dynamical attribute: they cluster 
just beyond the Hill stability boundary. Here we show that the planet-planet scattering model, 
which naturally explains the observed exoplanet eccentricity distribution, can reproduce the observed 
distribution of dynamical configurations. We calculated how each of our scattered systems would 
appear over an appropriate range of viewing geometries; as Hill stability is weakly dependent on 
the masses, the mass-inclination degeneracy does not significantly affect our results. We consider a 
wide range of initial planetary mass distributions and find that some are poor fits to the observed 
systems. In fact, many of our scattering experiments overproduce systems very close to the stability 
boundary. The distribution of dynamical configurations of two-planet systems actually may provide 
better discrimination between scattering models than the distribution of eccentricity. Our results 
imply that, at least in their inner regions which are weakly affected by gas or planetesimal disks, 
planetary systems should be "packed" , with no large gaps between planets. 
Subject headings: planetary systems: formation — methods: n-body simulations 



1. INTRODUCTION 

The observed eccentricities of extra-solar planets can 
be readily explained by a simple model that assumes 
that virtually all planetary systems undergo dynami- 
cal instabilities (Ford et al. 2003; Adams & Laughlin 
2003; Chatterjee et al. 2008; Juric & Tremaine 2008; 
Ford & Rasio 2008).* In the context of this model, 
planetary systems are expected to form in marginally 
stable configurations, meaning that they are stable for 
at least the timescale of rapid gas accretion of ^ 10^ 
years (Pollack et al. 1996) but ultimately unstable, prob- 
ably on a timescale comparable to the gaseous disk's 
hfetime of - 10^ years (Haisch et al.2001). This in- 
stability timescale implies an initial separation between 
planets of perhaps 4-5 mutual Hill radii Rh,m, where 
Rh,m ^0.5 {ai + a2)[{Mi+M2) /3M^]^/^; oi and 02 are 
the orbital distances. Mi and M2 are the masses of two 
adjacent planets, and is the stellar mass (Chambers 
et al. 1996; Marzari & Weidenschilhng 2002; Chatterjee 
et al. 2008).^ After a delay of 10*^ - lO'^ years, a typi- 
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* Several other models to explain the extra-solar eccentricity 
distribution exist; see Ford & Rasio (2008) for a summary. 

^ For Jupiter-mass planets, separations of ~ i — bRn m are close 
to the 3:2 and 2:1 mean motion resonances. Thus, an alternate 
argument in favor of planets forming with such spacings invokes 
resonant capture (Snellgrove et al. 2001) followed by turbulent re- 
moval from resonance (Adams et al. 2008) during the gaseous disk 
phase. 



cal system of three or more planets with separations of 
4 — 5Rh.m becomes unstable, leading to close encoun- 
ters between two planets, strong dynamical scattering, 
and eventual destruction of one or two planets by either 
collision with another planet, collision with the star, or, 
most probably, hyperbolic ejection from the system (Ra- 
sio & Ford 1996; Weidenschilhng & Marzari 1996; Lin & 
Ida 1997; Papaloizou & Terquem 2001). It is the plan- 
ets that survive the dynamical instability that provide a 
match to the observed extra-solar eccentricities. 

Additional dynamical information can be obtained 
from the known extra-solar multiple planet systems. In 
particular, the stability in two-planet systems can be 
guaranteed for planets with particular masses and orbital 
configurations. The edge of stability can be quantified in 
terms of the proximity to the analytically-derived Hill 
stability limit using the dimensionless quantity (3/Pcrit 
(the stability boundary is located at P/ficrit = 1; see 
Section 3). Dynamical analyses have shown that the 
known multiple- planet systems are clustered just beyond 
the edge of stabihty (i.e., at {3/ (icrit ^ 1; Barnes & Quinn 
2004; Barnes & Greenberg 2006, 2007). 

In this paper we study the stability of the surviving 
planets in several thousand 3-planet systems that have 
undergone planet-planet scattering leading to the loss of 
one planet. We find that in the aftermath of dynamical 
instabilities, the surviving planets cluster just beyond the 
stability boundary, providing a good match to the ob- 
served values. This provides support for planet-planet 
scattering as an active process in extra-solar planetary 
systems. This result also has consequences for the pack- 
ing of planetary systems and the "Packed Planetary Sys- 
tems" hypothesis (Barnes & Raymond 2004; Raymond 
& Barnes 2005; Raymond et al. 2006; Barnes et al. 2008). 
The paper proceeds as follows: we describe our scatter- 
ing simulations (§2), summarize Hill stability theory and 
define fi/Pcrit (§3), present our results (§4) and discuss 
the consequences (§5). 
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2. SCATTERING SIMULATIONS 

Our scattering simulations are drawn from the same 
sample as in Raymond et al. (2008a). Each simulation 
started with three planets randomly separated by 4-5 
mutual Hill radii. The three planets were placed such 
that the outermost planet was located two (linear) Hill 
radii Rh interior to 10 AU [Rh = a[M /iM^^/'^). We 
performed ten sets of simulations, varying the planetary 
mass distribution. For our two largest sets (1000 sim- 
ulations each) we randomly selected planet masses ac- 
cording to the observed distribution of exoplanet masses: 
dN/dM oc M^i i (Butler et al. 2006). In the Mixedl set 
we restricted the planet mass Mp to be between a Sat- 
urn mass Ms at and three Jupiter masses Mj^p- For our 
Mixed2 set, the minimum planet mass was decreased to 
10 M^. We also performed four Meq sets (500 simula- 
tions each) with equal mass planets for Mp = 30 M^, 
Msat, Mjup, and SMjup. Finally, the Mgrad sets (250 
simulations each) contained radial gradients in Mp. For 
the JSN set, in order of increasing orbital distance, Mp = 
Mjup, Msat, and 30 M®. For the NSJ set, these masses 
were reversed, i.e., the Mj^p planet was the most dis- 
tant. The 3JJS and SJ3J sets had, in increasing radial 
distance, Mp = 3Mj„p, Mjup and Msat, and Mp = Msat, 
Mjup and 3Mjup, respectively. 

Planetary orbits were given zero eccentricity and mu- 
tual inclinations of less than 1 degree. Each simulation 
was integrated for 100 Myr with the hybrid Mercury in- 
tegrator (Chambers 1999) using a 20 day timestep. We 
required that all simulations conserve energy to better 
than dE/E < lO"**, which is needed to accurately test 
for stability (Barnes & Quinn 2004). We achieved this 
by reducing the timestep to 5 days for simulations with 
dE/E > 10^^ and then removing simulations that still 
conserved energy poorly. As expected, these systems 
were typically unstable on 10^ — 10^ year timescales. In 
addition, about 1/4 of simulations were stable for 100 
Myr which shows that we started close to the stability 
boundary. For this paper, we restrict our analysis to the 
subsample of simulations that 1) were unstable, and 2) 
contained two planets on stable orbits after 100 Myr (i.e., 
one and only one planet was destroyed). 

3. HILL STABILITY 

For the case of two planets with masses Mi and M2 
orbiting a star, dynamical stability is guaranteed if: 

-2(Af^ + Mi + M2) 



1 + 34/3 



G^{MiM2 
M1M2 



M^Mi + M^M2) 
MiM2(llMi 



7M2 



iM^{Mi+M2Y 



■M2)4/3 

where c and h represent the total orbital angular mo- 
mentum and energy of the system, respectively (Marchal 
& Bozis 1982; Gladman 1993; Veras & Armitage 2004; 
note that this definition assumes that Mi > M2). We 
refer to the left side of Eqn 1 as /3 and the right side as 
Pcrit (Barnes & Greenberg 2006). The quantity (3/(3crit 
therefore measures the proximity of a pair of orbits to 
the Hill stability limit of P/Pcrit — 1- We note that our 
(3 / Peril analysis only applies for two-planet non-resonant 
systems, because perturbations from additional compan- 
ions can shift the stability boundary to values other than 
1 (Barnes & Greenberg 2007). 




Fig. 1. — Cumulative distribution of (i/f3crit of tlie well- 
ciiaracterized extra-solar multi-planet systems (in gray; see Table 
1), as compared with our scattering simulations. 

When calculating (3/ Pent for extra-solar systems, past 
research (Barnes & Greenberg 2006, 2007) has assumed 
coplanar orbits with masses equal to minimum masses. 
Those values of (3/ Pcrit were systematically affected by 
the mass-inclination degeneracy, probably resulting in 
overestimations. In contrast, our simulations provide the 
full three-dimensional orbits, and hence we can calcu- 
late the true value of P/Pcrit- More importantly, if we 
assume that viewing geometries are distributed isotropi- 
cally (i.e. edge-on systems are more likely than face-on), 
we can determine how P/ Pcrit would be calculated from 
radial velocity data (e.g. assuming coplanar, edge-on or- 
bits). For example, if two planets with masses and 
Mc have inclinations (relative to their invariable plane) 
ib and ic, and the inclination to the line of sight is /, then 
the "observed" P/ Pcrit would use masses Aff,sin(ib -I- /) 
and McSin{ic + /). In §4 we use this approach to build 
a distribution of P/ Pcrit that is directly comparable to 
the actual distribution (and effectively break the mass- 
inclination degeneracy) . 

4. RESULTS 

We generated P/ Pcrit distributions from our simula- 
tions following the procedure described above. First, we 
"observed" each system from 100 viewing angles, thereby 
decreasing the inferred mass of each planet by a factor 
of sin(/ + where ij refers to each planet's inclination 
with respect to a fiducial plane. Second, we assumed the 
observed systems to be coplanar in calculating P/ Pcrit 
for each viewing angle using Eqn 1 . Finally, we included 
the PI Pcrit calculated for each viewing angle by assum- 
ing the viewing angle / to be isotropically distributed. 
Figure [T] compares the cumulative P/ Pcrit distributions 
for the observed two-planet systems with our scatter- 
ing simulations. It is important to note that the "true" 
PI Pcrit distributions, calculated with knowledge of the 
simulated systems' real masses and inclinations, are vir- 
tually identical to the curves from Fig. [T] (this issue is dis- 
cussed further in §5). Table 1 lists the extra-solar systems 
in our analysis; we excluded systems with controversial 
or poorly-characterized orbits and those that were likely 
affected by tidal effects. In a two planet system with 
an inner planet at < 0.1 AU, tides will decrease the in- 
ner planet's eccentricity and semimajor axis (Jackson et 
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TABLE 1 

Planetary systems included in P/Pcru analysis-*- 









J vj ]_ • iv± 2j 




(pair) 


I" Ain 








HD 202206 b-c^ 


0.83,2.55 


0.435,0.267 


17.4,2.44 


0.883 


HD 82943 c-b^ 


0.746,1.19 


0.359,0.219 


2.01,1.75 


0.946 


HD 128311 b-c^ 


1.099,1.76 


0.25,0.17 


2.18,3.21 


0.968 


HI) 73^2fi h-r^ 


n 66 1 05 


0.19,0.14 


2.9 2.5 


0.982 




u.uo±,u.oy / 


n 1 fis n nQ7 
u. -i-Uo,u.uy 1 


n 1 S7 n fi^j^ 


u.yoy 




9 11 


n nAQ n 99 
u.u^y,u.zz 


9 ft n Afi 


1 09^ 


HD 155358 b-c 


0.628,1.224 


0.112,0.176 


0.89,0.504 


1.043 


HD 177830 c-b 


0.514,1.22 


0.40,0.041 


0.186,1.43 


1.046 


HD 60532 b-c2 


0.77,1.58 


0.278,0.038 


3.15,7.46 


1.054 


HD 183263 b-c 


1.52,4.25 


0.38,0.253 


3.69,3.82 


1.066 


HD 108874 b-c2 


1.051,2.68 


0.07,0.25 


1.36,1.018 


1.10 


HD 12661 b-c 


0.83,2.56 


0.35,0.2 


2.3,1.57 


1.12 


HD 11506 c-b 


0.639,2.43 


0.42,0.22 


0.82,3.44 


1.17 


HD 208487 b-c 


0.49,1.8 


0.32,0.19 


0.45,0.46 


1.20 


HD 169830 b-c 


0.81,3.60 


0.31,0.33 


2.88,4.04 


1.28 


HD 168443 b-c 


0.3,2.91 


0.529,0.212 


8.02,18.1 


1.95 


HD 38529 b-c 


0.129,3.68 


0.29,0.36 


0.78,12.7 


2.06 


HD 47186 b-c 


0.05,2.395 


0.038,0.249 


0.072,0.35 


6.13 



"See http: / / w-ww.astro.washington.edu / users /rory / research/xsp/dyi 
for an up to date list of f)/ [icrit values for the known extra-solar 
multiple planet systems. Orbital values were retrieved from 
|http://exoplanet.eu and http://exoplanets.org 

"These systems have been claimed to be in mean motion reso- 
nances. 

al. 2008), thereby increasing the separation between the 
two planets and l3//3crit- 

Four individual cases ~ Mixedl , Mixed2, Meq:Msat, 
and Meq:3OM0 - each provide a match to the observed 
/3/Pcrit distribution. Kolmogorov-Smirnov (K-S) tests 
show that the probability p that the (3/ Pcrit distributions 
from those four cases are drawn from the same distribu- 
tion as the observed sample are all 0.1 or larger (Table 
2). The distribution calculated by an unweighted com- 
bination of all ten cases is also a good match (each case 
was given equal weight, regardless of the number of sim- 
ulations). 

All of our sets of simulation produced a smaller fraction 
of systems at (3/ jScrit < 1 than for the observed systems. 
We therefore calculated K-S p values by confining the 
distributions to the ranges P/Pcrit < 1 and (3/ Pcrit > 1- 
All but one case with p > 0.1 also had p {P/Pcrit) > 0.1 
(Mixed2; see Table 2). However, some cases provide 
good matches for P/Pcrit < 1) but not for other regions, 
notably Meq:Mjup, Mgrad:JSN, and Mgrad:NSJ. Systems 
with P I Peril < 1 are unusual because they lie within the 
formal Hill stability boundary but are stabilized by spe- 
cial orbital configurations. In fact, all five of the known 
exoplanet systems with P / Pent < 1 are thought to lie in 
mean motion resonances (Table 1). The scattered sys- 
tems with P/Pcrit < 1 are stabilized by resonances or 
in many cases by low-amplitude, aligned apsidal libra- 
tion. Four cases in our sample generated resonant sys- 
tems in at least 5% of simulations (Raymond et al. 2008a) 
- Mixed2, Mgrad:JSN, Mgrad:NSJ, and Mgrad:SJ3J - 
but only two of these have p {P/Pcrit < 1) > 0.1. We 
attribute the lack of a correlation between resonances 
and P/Pcrit < 1 to the relative weakness of these res- 
onances. Indeed, resonances caused by scattering tend 
to exhibit relatively high-amplitude libration of only one 
resonant argument (Raymond et al. 2008a); these reso- 
nances have typical P/Pcrit values of slightly more than 



TABLE 2 

p values from K-S TESTS OF OBSERVATIONS VS. SCATTERING 
SIMULATIONS 





Case 


1 

P 




p{/3//3cr>t < 1) 


P {P/Pcr^t > 1) 


Mixedl 


0.14 




6.2 X 10-3 


0.53 


Mixed2 


0.13 




1.2 X 10""' 


0.07 


Meq:3Mjup 


1.8 X 10" 


6 


2.1 X 10-5 


4.6 X 10-" 


Meq:Mjup 


9.3 X 10" 


4 


0.10 


0.02 


Meq:Msat 


0.12 




0.81 


0.43 


Meq:3OM0 


0.29 




0.60 


0.89 


Mgrad: JSN 


4.0 X 10" 


3 


0.12 


2.4 X 10-" 


Mgrad:NSJ 


1.5 X 10" 


4 


0.10 


0.02 


Mgrad :3JJS 


9.1 X 10" 


3 


1.9 X 10-^ 


9.8 X 10-* 


Mgrad:SJ3J 


4.9 X 10" 


3 


6.2 X 10-4 


2.8 X 10-* 


All 10 cases 


0.14 




5.1 X 10-3 


0.58 



1 (median P/Pcrit — 1.01 — 1.03 for the different cases). 
This contrasts with resonances generated by convergent 
migration in gaseous disks, which tend to exhibit low 
amplitude libration of more than one resonant argument 
^SnpUgrove et al. 2001; Lee & Peale 2002). 

^Simulations with radial mass gradients (Mgrad) over- 
produced systems very close to the stability boundary, 
while cases with equal masses (Meq) produced much 
larger P/Pcrit values (Fig. [T|). A similar effect was seen 
in the eccentricity distributions: the Mgrad cases yielded 
much smaller eccentricities than the Meq cases (Raymond 
et al. 2008a; see also Ford et al. 2003). The Mixedl and 
Mixed2 cases fall between these two regimes. These 
trends can be explained by the number of close encoun- 
ters Uenc that occur in the different cases before a planet 
is destroyed. For the Mgrad typically be- 

tween 30 and 80, and is larger for less massive sys- 
tems (JSN and NSJ). For the Meq cases Ucnc is vastly 
larger, with median values between 100 (3Mjup) and 2000 
(30 Mq). The larger number of scattering events in- 
creases the eccentricity of surviving planets and also 
causes the systems to spread out farther. 

In calculating "observed" P/Pcrit distributions from 
our simulations, we assumed that the viewing angles I 
were isotropically distributed. Given that known extra- 
solar planet systems are each observed at a fixed I, could 
this have introduced a bias in our samples? Figure [2] 
shows the inferred value of P/Pcrit as a function of / 
for five Mixedl systems with varying mutual inclinations 
Ai. For Ai < 35°, P/Pcrit varies only slightly with the 
viewing angle, but for large Ai P/Pcrit can change sub- 
stantially with /.^"^ However, P/Pcrit varies by more than 
10% [20%] over the entire range of possible viewing an- 
gles for fewer than 10% [2%] of cases. For all systems, 
the edge-on P/Pcrit values agree with the true P/Pcrit 
values (calculated with knowledge of the planets' true 
masses and orbits) to better than 10%. There is a small 
bias: ~80% of systems exhibit a shallow negative slope 
in P/Pcrit vs. /, suggesting that the majority of inferred 
P/Pcrit values may be overestimated but only by < 1%. 
Thus, although / and Ai are important to keep in mind, 
they introduce a negligible error into our analysis. 

5. DISCUSSION 

We have found that it is actually the angular momentum 
deficit (Laskar 1997) which controls the magnitude of P/Pcrit vari- 
ation with /. 
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Ai = 




-^"^ Ai = 35° 




Ai = 21° - 




Al = 0° 




: 2:1 (9°) ; 



10 20 30 40 50 60 70 80 90 

I (°) 

Fig. 2. — Inferred values for /3//3crit as a function of observation 
angle / for several examples from the Mlxedl set, labeled by the 
approximate mutual inclination Ai between planets. One resonant 
case is labeled "2:1": it has Ai = 9°. / = 90° is edge-on and 
7 = 0° is face-on. 

The planet-planet scattering model appears to be con- 
sistent with the {3/ Peril distribution of the observed ex- 
oplanct systems. The distribution can be reasonably re- 
produced by several of our sets of simulations, or even by 
an unweighted combination of all ten sets. We therefore 
cannot strongly constrain the initial planetary mass dis- 
tribution, although we can rule out cases with very poor 
fits Meq:3Mjup, Mgrad:3JJS and Mgrad:SJ3J in partic- 
ular as the major contributors to the distribution (see 
Table 2). We consider the Mixedl set to be the most re- 
alistic because it is drawn from the observed mass distri- 
bution (Butler et al. 2006), and it provides a good match 
to the observed eccentricity distribution (Raymond ct 
al. 2008a). In the coming years, we expect many more 
systems to be discovered with P/ ficrit « 1 — 1-5. 

The pileup of scattered systems just beyond the 
stability boundary implies that planetary systems are 
"packed", meaning that large spaces in between plan- 
ets should be rare.^^ This provides a theoretical foun- 
dation for the "Packed Planetary Systems" hypothesis, 
which asserts that if a stable zone exists between two 
known planets, then that zone is likely to contain a planet 
(Barnes & Raymond 2004; Raymond & Barnes 2005; 
Raymond et al. 2006). Given the small /3//3crit values 
of scattered systems, there is simply no room to insert 

It is important to note that the spacing for planets with 
PI Pcrit ^ 1 can be large. Among just the Mixedl simulations with 
1 ^ P/Pcrit < 1-1 the difference in semimajor axis for adjacent 



another planet between the two known planets without 
causing the system to be unstable. 

HD 74156 is an example of a packed planetary sys- 
tem. Prior to 2008, two planets were known in the 
system, at 0.28 and 3.4 AU (Naef et al.2004), with 
(i/Pcrit = 1-987. Raymond & Barnes (2005) mapped 
out a narrow stable zone between the two planets, from 
0.9-1.4 AU. The planet HD 74156 d was discovered three 
years later by Bean et al. (2008) at 1.01 AU (see also 
Barnes et al. 2008) at the peak of the stable zone. We 
therefore expect additional planets to exist in systems 
with PI Pent > 1.5-2, notably HD 38529 (Raymond 
& Barnes 2005) and HD 47186 (Kopparapu ct al. 2009). 
The probable location of additional planets can be deter- 
mined using test planets to map out dynamically stable 
regions between known planets (e.g., Menou & Tabach- 
nik 2003: Rivera & Haghighipour 2007; Raymond et 
al. 2008bj. 

The P/Pcrit distribution of the observed extra-solar 
planetary systems may contain information about dif- 
ferent dynamical regimes. The region of P/Pcrit < 1 
is populated entirely by resonant systems and may pro- 
vide evidence of planetary system compression, presum- 
ably via convergent migration in gaseous protoplane- 
tary disks (Snellgrove et al.2001; Lee k Peale 2002). 
The region of 1 < P/Pcrit < 1.5 — 2 is consistent with 
the scattering regime. Widely-separated systems with 
P/Pcrit > 1.5 — 2 may have been drawn apart by inter- 
actions with planetcsimal or gaseous disks (e.g.. Gomes 
et al.2004; Moeckel et al. 2008). However, this seems 
unlikely given that the known planets lie relatively close 
to their stars and that disk effects should be far more 
pronounced at large distances. 

Given that our simulations started with only three 
planets, we could not calculate P/Pcrit values in per- 
turbed two-planet systems. For example, an interest- 
ing comparison with observations would be to measure 
P/Pcrit for the two easiest-to-detect planets in scattered 
three planet systems. This would address the question 
of whether to search for additional planets in between 
or interior/exterior to the known planets in two planet 
systems with large P/Pcrit- 

We thank Google for access to their machines. S.N.R. 
and R.B. acknowledge funding from NASA Astrobiol- 
ogy Institutes 's Virtual Planetary Laboratory lead team, 
supported by NASA under Cooperative Agreement No. 
NNH05ZDA001C. 



planets ranges from <2 AU to >15 AU. 
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